Sensor vs Monitor

een analyse van drie locaties voor PM_2,5-metingen in de stad Antwerpen voor de periode 2021-22.

Inleiding

In dit document doen we een korte vergelijkende studie tussen de fijn-stof-metingen van drie van de VMM meetstations in de stad Antwerpen met de low-cost PM_2,5 sensoren die op dezelfde locaties geinstalleerd staan. We zullen ons focussen op drie meetplaatsen met locaties als volgt:

Deze kaart toont de locatie van de drie meetstations die we in deze analyse in rekening nemen. In wat volgt zullen we naar hen verwijzen met hun straatnaam; zijnde de Plantin en Moretuslei (P), de Belgiëlei (B) en de Groenenborgerlaan (G). We tonen bovendien ook de locatie van het KMI weerstation te Stabroek waarvan we data gebruiken om extra context te verschaffen over onze metingen.

In een latere sectie van dit document zullen we tevens weerdata betrekken bij onze analyse, komende van het geautomatiseerd weerstation van de KMI te Stabroek. Met behulp van deze informatie hopen we de discrepantie tussen sensor- en monitordata te contextualiseren.

Sensordata

We beginnen door de sensordata voor de verschillende locaties samen te plotten voor de jaren 2021-2022. Zo krijgen we een eerste indruk van de structuur van de meetresultaten, alsook kunnen we een aantal kwantitatieve kenmerken vaststellen. De onderstaande plot is een interactieve, waarop ingezoomd kan worden naar een specifieke periode naar keuze. Het is van belang om te vermelden dat we verscheidene metingen voor de meetplaats op de Plantin en Moretuslei verwijderd hebben omdat ze boven de \(500 \, \frac{\mu g}{m^3}\) lagen en dus weggezet kunnen worden als anomaliën.

Een eerste kwestie die onmiddelijk opvalt, is dat de tijdsreeksen voor de sensoren op regelmatige basis onderbroken worden. Vermoedelijk is dit te wijten aan verbindings-problemen of andere defecten en zijn dit onvermijdelijke aspecten van het gebruiken van low-cost automatische metingen. Om de mancos in de data te visualiseren, plotten we hieronder de tijdsreeksen apart met zwarte banden voor de periodes zonder metingen.

In totaal ontbreken er voor de meetpunten in de Plantin en Moretuslei 9%, voor de Belgiëlei 6% en voor de Groenenborgerlaan 8% van de datapunten.

Een tweede feit dat we onmiddelijk kunnen opmerken, is dat de data vrij erratisch verloopt en scherpe pieken vertoont. Als we de daggemiddelden plotten per sensor krijgen we een iets grover beeld dan de hierboven getoonde uurlijkse metingen.

Het mag weinig verbazen dat de drastische pieken hiermee niet volledig uitgemiddeld worden, maar dit geeft ons wel een overzichtelijker beeld van hoe de locaties zich onderling verhouden. We zien dat ze over het algemeen vrij gelijk verlopen, maar dat er occasionele pieken zijn voor individuele sensoren die niet te zien zijn in de andere metingen. Vergelijken met de monitordata zou ons moeten kunnen leren of dit te wijten is aan hyperlokale effecten, dan wel aan tijdelijke storingen aan het meetinstrument.

Laten we dan, om robuustere grootheden dan de daggemiddelden te bestuderen, een gemiddeld dag-, week- en maandverloop bekijken. Hiervoor kunnen we standaard-functies uit het R-pakket openair gebruiken. Deze plots tonen een samenvatting van de gemiddelde bovengenoemde periodes, samen met een 95% getrouwheidsinterval. Zo krijgen we ook een beeld van de spreiding van de data.

Wat door deze plots onmiddelijk in het oog springt, is dat de gemiddelde week- en dagverlopen relatief gelijk lopen tussen de verschillende sensoren. Dit ligt in de lijn der verwachtingen, aangezien de meetlocaties dicht bij elkaar liggen, en we zouden verwachten dat een grootheid die robuuster is voor erratische schommelingen ook gelijkaardig zou moeten zijn voor de verschillende sensoren, aangezien de relevante afhankelijken - verontreinigingsbronnen en meteorologische data - relatief gelijkaardig zouden moeten zijn.

Verschil met Monitordata

Laten we dezelfde plots genereren voor de monitordata in de hoop hierdoor al enige verschillen te kunnen waarnemen.

Het valt dadelijk op dat de monitordata veel sterker gelijk loopt tussen de verschillende locaties. We durven dus vermoeden dat de occasionele grote verschillen tussen locaties onderling voor de sensoren het gevolg waren van meetfouten.

Als we de daggemiddelen plotten is deze strakke overeenkomst nog vele malen duidelijker. Reeds deze eerste blik levert dus al een onderscheid tussen monitor en sensordata op. Om dit expliciet te maken, definieren we voor elke locatie een dataframe dat zowel de sensor- als de monitordata bevat, en berekenen we hun verschil via \[\rho^{PM_{2.5}}_{diff} = \rho^{PM_{2.5}}_{sensor} - \rho^{PM_{2.5}}_{monitor}.\] Laten we de tijdsreeks voor de verschillen voor elke locatie plotten

en tevens de volgende beschrijvende getallen berekenen:

• Gemiddelde Absolute Fout (MAE): Meet het gemiddelde van de absolute verschillen tussen sensor- en monitorgegevens.
• Gemiddelde Kwadratische Fout (MSE): Berekent het gemiddelde van de gekwadrateerde verschillen tussen sensor- en monitorgegevens.
• Wortel van de Gemiddelde Kwadratische Fout (RMSE): RMSE is de vierkantswortel van de MSE.
• Gemiddelde Afwijkingsfout (MBE): Meet het gemiddelde verschil tussen sensor- en monitorgegevens.
• Standaarddeviatie van de Verschillen (SD): De standaarddeviatie van de verschillen.
• Mediaan Absolute Afwijking (MAD): meet de mediaan van de absolute verschillen tussen sensor- en monitorgegevens.
• Maximale en Minimale Verschillen (Max/Min): de maximale en minimale verschillen.
• Interkwartielafstand (IQR) van Verschillen: vertegenwoordigt het bereik waarbinnen de middelste 50% van de verschillen ligt.
• Variantiecoëfficiënt (CV) van Verschillen: drukt de relatieve variabiliteit van de verschillen ten opzichte van het gemiddelde verschil uit.

Table 1: Beschrijving Verschil Sensor en Monitor

Locatie	MAE	MSE	RMSE	MBE	SD	MAD	Max	Min	IQR	CV
Plantin en Moretuslei	6.09	73.90	8.60	-3.47	7.87	3.31	98.67	-205.89	4.47	-2.27
Groenenborgerlaan	5.55	81.51	9.03	-0.43	9.02	3.76	124.36	-151.96	5.58	-20.94
Belgiëlei	5.27	54.06	7.35	-2.63	6.87	3.16	109.74	-191.72	4.38	-2.62

De negatieve MBE wijst erop dat de sensoren overwegend lijken te onderschatten. Bovendien merken we disproportioneel hoge waarden voor de minimum en maximum verschillen, desondanks valt de omvang van de verschillen voor \(50%\) van de tijd binnen het interval \(6.09 \pm 4.47\), \(5.55 \pm 5.58\) en \(5.27 \pm 4.38\) voor de drie locaties. Om het onderscheid tussen de monitor- en sensordata verder duidelijk te maken, genereren we een scatterplot met lineaire regressie voor elke locatie.

Table 2: Lineaire Regressie-coefficienten

Locatie	Intercept	Slope	R²
Plantin en Moretuslei	-0.5276234	0.7661490	0.5022223
Groenenborgerlaan	-0.1229899	0.9730239	0.4956407
Belgiëlei	-1.0566014	0.8822739	0.6424049

In het ideale scenario zouden de monitor en de sensor in alle gevallen moeten overeenstemmen en dus een een-op-een relatie hebben. Aan de scatterplots zien we echter duidelijk dat dit niet het geval is. Wanneer we de coefficienten van de lineaire regressie overlopen zien we niet enkel afwijkingen in het intercept en de regressiecoefficient, maar vooral een beduidend lage determinatiecoefficient. In deze zien we dus (zoals reeds visueel duidelijk van de scatter plot) dat monitor en sensor geen eenduidig lineair verband hebben.

Wat zijn de storende factoren die dit verband bemoeilijken? Meteorologische factoren zijn van belang voor de accuratie van low-cost PM sensoren #TODO: goeie citatie vinden Het is bijvoorbeeld geweten dat vele van de componenten van fijn stof hygroscopisch zijn. Hierdoor speelt bijvoorbeeld luchtvochtigheid een rol in het interpreteren van de sensormetingen, aangezien het water dat de stofdeeltjes aantrekken de meting zelf kan gaan verstoren. Vele studies hebben reeds aangetoond dat het in rekening brengen van meteorologische factoren bij het verder kalibreren van sensordata tot goede resultaten kan leiden.

Laten we dus data van het KMI weerstation voor de relevante periodes toevoegen aan onze dataset in de hoop hiermee het verschil tussen de sensor- en monitordata verder te gaan interpreteren en (mogelijks) verkleinen. We voegen aan onze dataframes relevante meteorologische informatie zoals luchtvochtigheid (RH), neerslag (RR), windrichting (wd) en windsnelheid (ws) toe. Om hun mogelijke impact op \(\rho^{PM_{2.5}}_{diff}\) vast te stellen, berekenen we de correlatiematrices voor deze variabelen.

We zien dat deze resultaten het vermoeden bevestigen dat meteorologische grootheden een impact hebben op het verschil tussen sensor- en monitordata. In het bijzonder de luchtvochtigheid, maar ook windsnelheid (voor de Groenenborgerlaan) lijken een rol te spelen. Om deze afhankelijkheden mede in rekening te nemen bij het analyseren van de sensordata, volgen we het voorbeeld van e.g. (Bush et al. 2022) en passen we in de volgende sectie machine learning technieken toe om deze effecten deels te corrigeren.

Kalibratrie van Sensordata via RF Regressie

In navolging van verscheidene papers uit de literatuur (deSouza et al. 2022; Bush et al. 2022; Wang et al. 2019), gebruiken we een Random Forest (RF) regressiemodel, een methode die efficient en robuust blijkt (Liang and Daniels 2022). Dit is een zogenaamd ensemble learning model dat gebruik maakt van een grote verzameling beslisbomen, elk opgesteld voor een willekeurige deelverzameling van de oorspronkelijke test-data. Door uit te middelen over een groot aantal beslisbomen is deze methode robuuster tegen overfitting dan de traditionele beslisboom-methodes.

In ons geval gebruiken we de volledige sensor-tijdsreeks van 2021 en 2022 als testdata, en passen we een RF model toe om een regressie hierop uit te voeren voor elke locatie, om zo een post-hoc gekalibreerde sensor te verkrijgen die hopelijk dichter bij de monitordata aanleunt. De details van de modellen (beslissende parameters zoals aantal bomen) werden getweaked om een optimaal resultaat te bekomen en vindt men hieronder opgelijst

Call:
 randomForest(formula = pm25 ~ ., data = bigframe_p, ntree = 2000,      mtry = 4, importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 2000
No. of variables tried at each split: 4

          Mean of squared residuals: 16.10012
                    % Var explained: 83.46

Call:
 randomForest(formula = pm25 ~ ., data = bigframe_g, ntree = 2000,      mtry = 4, importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 2000
No. of variables tried at each split: 4

          Mean of squared residuals: 17.61472
                    % Var explained: 79.14

Call:
 randomForest(formula = pm25 ~ ., data = bigframe_b, ntree = 2000,      mtry = 4, importance = TRUE) 
               Type of random forest: regression
                     Number of trees: 2000
No. of variables tried at each split: 4

          Mean of squared residuals: 18.37986
                    % Var explained: 82.59

De modellen die hierboven staan opgelijst genereren dus een gekalibreerde versie van de sensordata die poogt de monitordata te benaderen. We kunnen net als voorheen een scatter plot maken die deze onderling vergelijkt, en een lineaire regressie uitvoeren.

Table 3: Lineaire Regressie-coefficienten Gekalibreerde Sensor

Locatie	Intercept	Slope	R²
Plantin en Moretuslei	2.548079	0.8004830	0.8361167
Groenenborgerlaan	2.653228	0.7723842	0.7919514
Belgiëlei	2.593157	0.8088849	0.8263124

Hiermee zien we onmiddelijk dat de gekalibreerde sensordata voor 2021 en 2022 aanzienlijk dichter aanleunt bij de monitordata. Dit ligt echter in de lijn der verwachtingen, aangezien we deze periode ook gebruikten als testdata. Hoewel het intercept eigenlijk vergrootte en de regressiecoefficienten nagenoeg constant bleven, zien we vooral een aanzienlijk hogere determinatiecoefficient, en is de lineaire regressie dus aanzienlijk betrouwbaarder. Dit impliceert dat de kalibratie een groot deel van de ongewenste storingsfactoren kan wegnemen.

We kunnen wederom het verschil tussen de gekalibreerde sensor en de monitor gaan berekenen voor elke meetgebeurtenis. Indien we hiervan enkele beschrijvende cijfers berekenen voor elke locatie, krijgen we de volgende waarden.

Table 4: Beschrijving Verschil Gekalibreerde Sensor en Monitor

Locatie	MAE	RMSE	MBE	SD	MAD	Min	Max	IQR	CV
Plantin en Moretuslei	2.30	4.01	0.04	4.01	1.97	-156.43	52.40	2.66	100.25
Groenenborgerlaan	2.31	4.20	0.05	4.20	2.00	-131.14	69.27	2.69	84.00
Belgiëlei	2.33	4.29	0.04	4.29	2.01	-170.09	142.05	2.72	107.25

Hiermee wordt dan onmiddelijk duidelijk dat de kloof tussen de sensor en monitordata dankzij de RF regressie aanzienlijk gekrompen is. Dit valt binnen de lijn der verwachtingen voor de testdata, laten we nu de data voor 2023 ophalen om te zien of ons model buiten de test-range nuttige resultaten oplevert.

RF Model toepassen op 2023 Data

We importeren de sensor- en monitordata van 01-01-2023 tot 06-09-2023 en passen hierop dezelfde regressie toe. Eerst en vooral vergelijken we zoals voorheen de ongekalibreerde sensor met de monitordata en zien in hoeverre hun onderling verband afwijkt van een simpele lineaire kromme.

Table 5: Lineare Regressie-coefficiënten Sensor 2023

Locatie	Intercept	Slope	R²
Plantin en Moretuslei	-2.3582003	0.8380143	0.7323903
Groenenborgerlaan	-2.0401687	1.2476051	0.4097593
Belgiëlei	-0.6906443	0.7214090	0.6306911

We zien wederom dat er enige ruis op de overeenkomst tussen monitor en sensor zit. Om een beter idee te krijgen van hun verschil geven we wederom een tabel van enkele beschrijvende waarden voor \(\rho^{PM_{2.5}}_{diff}\) in de relevante periode.

Table 6: Beschrijving Verschil Sensor en Monitor 2023

Locatie	MAE	RMSE	MBE	SD	MAD	Min	Max	IQR	CV
Plantin en Moretuslei	5.24	6.51	-4.29	4.90	3.19	-28.28	35.10	4.36	-1.14
Groenenborgerlaan	5.23	10.49	0.36	10.48	3.35	-19.35	140.71	4.74	28.89
Belgiëlei	4.98	6.97	-4.04	5.69	3.72	-165.25	35.96	5.07	-1.41

Ook hier zien we dat de verschillen tussen de sensor en monitordata gelijkaardig zijn in 2023. Een interessant detail we voor het eerst een positieve MBE waarnemen op de Groenenborgerlaan, wat impliceert dat de sensor hier systematisch overschat. Het valt te bezien of de oorsprong hiervoor ligt in problematisch gedrag van de niet-gevalideerde monitordata.

Nu passen we de RF modellen toe op de 2023 data, en vergelijken we de gekalibreerde sensordata met de monitordata.

Table 7: Lineaire Regressie-coefficienten Gekalibreerde Sensor 2023

Locatie	Intercept	Slope	R²
Plantin en Moretuslei	4.150791	0.6727324	0.7782443
Groenenborgerlaan	3.275557	0.9083790	0.5569169
Belgiëlei	4.385430	0.6310281	0.6109426

Table 8: Beschrijving Verschil Gekalibreerde Sensor en Monitor 2023

Locatie	MAE	RMSE	MBE	SD	MAD	Min	Max	IQR	CV
Plantin en Moretuslei	2.79	4.48	0.26	4.48	2.44	-40.91	13.78	3.30	17.56
Groenenborgerlaan	3.36	6.11	2.39	5.63	2.22	-30.03	105.46	3.03	2.36
Belgiëlei	3.84	5.74	-0.05	5.74	4.18	-158.41	17.44	5.64	-124.17

We zien dat onze RF regressie weliswaar een positieve impact heeft op de omvang van de discrepantie tussen monitor en sensor voor de 2023 data, maar dat de lineaire regressie er niet op verbetert na kalibratie. Dit doet vermoeden dat ons model nog niet optimaal is en er werk aan de winkel is om deze resultaten te verbeteren.

Enkele mogelijke volgende pistes:

1. De monitordata die we momenteel gebruiken voor 2023 is niet gevalideerd, het zou dus kunnen dat er problematische datapunten in voorkomen die onze analyse verstoren.

2. We hebben onze tijdsreeks niet gedecorreleerd vooraleer het model te trainen, mogelijks is dit een piste die bekeken moet worden.

3. In verscheidene papers gebeurt er eerst een baseline correction vooraleer over te gaan op het regressiemodel. Mogelijks kan dit de efficientie verhogen.

4. In analyses zoals (Adong et al. 2022) werd ook temperatuur toegevoegd aan het regressiemodel, mogelijks missen we dus nog deze variabele.